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ABSTRACT 


A  comparison  of  several  techniques  is  presented  for 
determining  upper  confidence  levels  for  a  system  failure 
rate.   A  series  system  of  components  with  exponential 
failure  rates  is  examined.   Classical  computational  tech- 
niques are  compared  with  Bayesian  techniques  in  determining 
the  upper  confidence  level  of  a  system  failure  rate.   A 
sensitivity  analysis  is  conducted  on  several  of  the 
parameters  as  part  of  the  comparison. 
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I.   INTRODUCTION 

A.   BACKGROUND 

Numerous  "classic"  techniques  have  been  used  to  compute 
estimates  of  failure  rate  and  mean  time  of  failure.   From 
these  estimates  standard  accepted  procedures  can  be  applied 
to  establish  upper  confidence  levels  (UCL)  for  the  failure 
rate  or  lower  confidence  levels  (LCL)  for  the  reliability. 

One  well  established  "classic"  procedure  is  to  utilize 
the  computational  methods  set  forth  in  Ref.  3  to  determine 
upper  confidence  levels  on  system  failure  rates.   This  is 
the  procedure  used  for  all  "Classic"  and  "Semi-Classic" 
methods  presented  in  this  paper. 

The  application  of  a  "Bayesian"  technique  may  be  intui- 
tively appealing  to  some.   Results  from  previous  experiments 
and  testing  could  be  applied  apriori  to  current  testing 
programs  to  determine  failure  rate  and  reliability.   By- 
using  a  prior  in  such  a  manner  perhaps  total  system  test 
time,  number  of  component  failures,  etc.,  could  be  reduced, 
thereby  reducing  the  overall  expense  involved  with  system 
life  testing.   This  would  appear  to  be  most  appealing  when 
testing  systems  that  are  extremely  expensive. 

The  intent  of  this  thesis  is  to  attempt  to  determine  if 
it  would  be  more  advantageous  to  use  some  "Bayesian"  tech- 
nique rather  than  a  more  traditional  "classic"  technique 
to  compute  upper  confidence  levels  on  a  system  failure  rate 


when  the  "prior"  for  the  Bayesian  method  could  be  chosen 
as  "optimistically"  as  one  would  desire. 

B.   SYSTEM  ASSUMPTIONS 

Since  it  is  the  technique  of  computation  of  the  upper 
confidence  level  that  is  being  investigated,  the  system 
that  is  being  modeled  can  be  kept  simple,  yet  still 
realistic. 

Even  the  most  complicated  system  can  be  broken  down  to 
a  system  of  components,  connected  in  series,  which  the 
total  failure  of  any  one  component  will  result  in  the  mis- 
sion failure  of  the  system. 

Each  type  component  in  this  series  system  experiences 
exponential  failure  rate.   The  failure  of  each  type  of 
component  is  assumed  to  be  independent  of  the  failure  of 
any  other  type  of  component.   Wearin  and  wearout  are  neg- 
lected and  the  failure  rate  is  assumed  to  be  constant. 
Components  are  assumed  not  to  ever  be  "stillborn."  To 
keep  the  system  as  simple  as  possible,  each  type  component 

exhibits  an  identical  failure  rate,  X,  -  X.;  i  =  2  .  .  .  k, 

*   1    1 ' 

k  equaling  the  number  of  type  components  that  are  connected 
in  the  series  system.   The  following  formulas  concerning 
system  failure  rate  and  system  reliability  are  assumed  to 
be  valid  for  this  system: 

k 
A   =  Z   A- 
5    i  =  l  x 


R   =  R-l     or  R  =  exp~Xs 

In  keeping  the  system  as  simple  as  possible  the 
optimistic  "prior"  for  the  Bayesian  computations  are  also 
chosen  to  be  equal  for  each  type  component,  a-  «  a.  and 
3   =  3- ,  for  i  =  1,  2,  .  .  .  k. 

For  purposes  of  comparison,  it  is  assumed  that  each 
type  of  component  is  placed  on  test  for  an  equal  length  of 
time,  t,  =  t.,  and  each  type  of  component  experiences  the 
same  number  of  failures  during  its  total  test  time, 
f,  =  f • .   This  assumption  is  modified  slightly  to  be  able 
to  examine  the  system  that  has  only  one  component  failure. 
A  system  exhibiting  two  component  failures  is  also 
examined. 

The  data  (values  of  the  parameters)  have  been  chosen 
by  the  author.   They  have  been  intentionally  chosen  to  be 
computationally  simple  yet  to  still  exhibit  the  characteris 
tics  of  a  realistic  system. 


II.   COMPUTATIONAL  METHODS 

A.   CLASSIC 

Reference  3  is  used  exclusively  for  the  computations  in 
this  method.   The  computational  formulas  given  in  this 
reference  are  modified  slightly  to  accommodate  the  basic 
assumptions  of  identical  test  time,  identical  number  of 
failures  for  each  component,  etc. 

X.  is  a  Maximum  Likelihood  Estimate  for  the  failure 

1 

rate  of  the  i   component.   It  will  be  equal  to  the  total 
number  of  failures  divided  by  the  total  test  time. 

X  ,  the  system  failure  rate,  will  be  equal  to  the  sum 

of  the  individual  component  failure  rates.   X   =   EX.. 

s    i  =  l  1 

Each  component  is  assumed  to  fail  independently  and  k  w:  .1 

equal  the  number  of  components  that  are  placed  together 

in  a  series. 


X  = 
u 


?       2       ~    2       421/2 
2Xc  +  K'^C  +  (4XK,ZC  +  K'4C^)1/Z 


X  will  be  an  upper  confidence  level  on  the  system 

failure  rate. 

k  X. 

c  =  izl-i 


The  iralues  for  K'  are  found  in  Table  I.   The  computational 
procedure  to  determine  the  Beta  values  for  a  90  per  cent 
level  of  significance  are  also  discussed.   Eighty  per  cent 
level  of  significance  values  are  found  in  Ref.  3. 

A  formula  for  the  upper  confidence  level  for  the  system 
failure  rate  is  also  provided  when  no  failures  have  occurred 
during  the  total  test  time  that  each  component  has  been 
allotted.   That  formula  follows: 


,2   k 
Z 
i  =  l   "i 


ILL.    '  1 

u    n 


i   (^ 


where  n  equals  the  number  of  component  terms  in  ths 
summation. 

Substitution  of  the  upper  limits  on  the  failure  rates 
obtained  will  generate  corresponding  lower  confidence 
limits  on  reliability*,  thus 


Rn  =  exp"Xu  . 


Beta  can  be  calculated  from  the  following  formula 


X   -  X  =  Beta(K) (X  )1/2 
u  v  J  ^   uJ 


Xu  is  obtained  from  Chi-Squared  tables  in  Ref.  1  at 
the  desired  confidence  level  with  2(X+1)  degrees  of  free 
dom  where  X  is  the  number  of  failures  and  K  is  the 
percentage  point  of  the  normal  distribution  point  in  the 
same  confidence  level. 
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TABLE  I 
BETA  VALUES  FOR  90  PER  CENT  CONFIDENCE 

Mr.   r.-F    y2  m/?y2    v  ~\      Beta        K1 

ino   or    a   90  2(X+1)   ^1/ZA  .90,2(X+1)J   Value   (Beta  (1 .  282) ) 
Failures       '  *■    '  '    v    '  K  K  JJ 

0  4.906  2.3025       1.18362    1.5174 

1  7.779  3.889        1.14271    1.46496 

2  10.645  5.3225  1.12336  1.44015 

3  13.362  6.681  1.1109  1.42417 

4  15.987  7.994  1.10188  1.41261 

5  18.549             9.275  1.09494  1.403713 
10  30.813  15.406  1.0743  1.37725 
20  54.090  27.045  1.05669  1.35468 
30            INTERPO  L.A  TED  1.34639 
36.5  91.1  45.55  1.04596  1.34092 

The  "Beta  Value"  listed  in  this  table  is  germane  to 
Ref.  3  and  should  not  be  confused  with  the  3  that  is  the 
optimistic  scale  parameter  used  in  the  Bayesian  simulation. 

B.   BAYESIAN 

The  Bayesian  technique  will  assume  the  apriori  density 

to  be  GammafX . ;a. ,3 • ) •   a •  and  3-  are  the  shape  and  scale 
v  i   i   l     i      l  r 

parameters.  The  "prior"  can  be  made  "optimistic"  if  the 
scale  parameter  is  chosen  large  in  relation  to  the  shape 
parameter. 
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The  a  posteriori  density  then  becomes  (A . ;a. +£.,$. +t .) . 

f.  will  equal  the  number  of  failures  of  the  ith  type  com- 

ponent  and  t.  will  equal  the  total  test  time  for  the  i 

type  component. 

k 
The  distribution  of  X   =   EX.  will  then  be  determined 

s    i  =  l  x 

by  computer  simulation.   Random  variates  of  each  X.  are 

generated  from  the  Gamma  distribution  with  parameters 

a.  +  f - ,  3-  +  t . .   A  random  variate  will  be  generated  for 
1111  to 

each  X.,  i  =  1,  .  .  .  k,  the  number  of  components  in  series. 

The  X..   will  then  be  added  to  determine  the  series  system 

i's  ' 

failure  rate. 

The  process  for  generating  a  system  failure  rate  is 

then  repeated  1000  times  to  yield  X  ,,  X  „,  .  .  .  X  nnnn. 
r  3  si'   s2'        slOOO 

The  1000  random  values  of  X   are  then  ordered  to  yield 

s 

Xs(l)'  Xs(2)'  *  *  *  Xs(1000)* 

■f"  v> 
X   (  ->  is  the  "estimate"  of  the  (1  -  y)    percentile 

point  of  the  distribution  of  X  . 
r  s 

The  Bayesian  100(1  -  y)    upper  confidence  level  for 

X„  is  then  X  -.nnnr-i  ^ 

s  sl000(l  -  yj 

At  the  time  of  this  writing  a  subroutine  to  generate 

random  variates  from  the  Gamma  distribution  does  not  exist 

in  the  Naval  Postgraduate  School  computer  library.   Reference 

2  was  used  to  generate  random  variates  from  the  Gamma 

distribution  with  an  integer  shape  parameter. 
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The,  generation  of  random  variat.es  with  a  shape  parameter 
that  is  not  an  integer  poses  a  much,   more  complicated  prob- 
lem.  Reference  4  has  been  written  to  handle  this  situation 
for  shape  parameters  between  0..05  and  1.0.   Although 
LT  Robinson's  approach  is  relatively  untried,  it  shows  a 
great  deal  of  promise  and  may  possibly  be  incorporated 
into  the  Naval  Postgraduate  School  computer  library  in  the 
future.   It  has  been  used  in  this  writing  for  comparative 
purposes  for  a  shape  'parameter  less  than  1.0. 

C.   SEMI-CLASSIC 

Similar  to  the  Classic  technique,  Ref.  3  is  used  to 
compute  the  upper  confidence  levels  for  the  system  failure 
rate  for  this  method.   The  computations  will  be  modified 
somewhat  by  adding  the  identical  "optimistic"  priors  used 
in  the  Bayesian  technique. 

The  number  of  failures  per  component  will  be  added  to 
a.  and  the  total  test  time  per  component  will  be  added  to 
3-.   A  maximum  Likelihood  Estimate  for  the  i    component's 
failure  is  then: 

f.  +  a. 

The  system  failure  rate  for  this  Semi-Classic  method  is 
the  same  as  for  the  Classic  method. 

xs  =  *h 
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The  value  of  C  will  also  change  slightly 

.A, 

x. 


t   +8 
C  =   1    pi 


Xs 

When  no  failures  have  occurred,  the  upper  confidence  limit 
on  the  failure  rate  now  is: 


2 
X   =  K' 


u   n     t.  +  8- 
i    l 

The  computation  for  the  upper  confidence  level  with  failure 
remains  the  same. 

This  method  is  somewhat  "Bayesian"  in  the  sense  that 
it  is  utilizing  a  prior  but  it  is  "Classic"  because  of  the 
nature  of  the  computational  method; 


14 


III.   PARAMETERS  OF  VARIATION 

Although  the  range  through  which  some  of  the  parameters 
are  varied  could  be  more  extensive,  it  does  provide  the 
reader  with  an  adequate  base  for  comparing  the  methods  of 
computation. 

a.       1.0,0.5 

3.      50,  000 

i        ' 

t.       30,  50,  100  mission  units 

f.      1.0,  0.0   The  failures  are  also  modified  so 

that  f,  =1,  f~  =  f.  =0  and  f,  =  1,  f0  =  1, 
1     *   2    i  1     '   2    ' 

f.,  =  f .  ,  f.'  =  0.   This  will  provide  the 

reader  with  the  opportunity  to  compare  a 
system  where  only  one  component  or  only  two 
components  fail. 

k       1,  2,  5,  10,  30  components 

Y       0.10 
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IV.   COMPARISONS  OF  SYSTEM  FAILURE  RATE 
UPPER  CONFIDENCE  LEVELS 

The  following  12  tables  provide  the  reader  with  upper 
confidence  levels  computed  from  the  same  data  by  the  three 
different  methods:   Classic,  Semi-Classic,  and  Bayesian. 

A  program  to  generate  random  variates  from  a  Gamma 
distribution  with  a  shape  parameter  that  was  greater  than 
1.0  but  not  an  integer  was  not  available;  therefore  UCL's 
for  the  Bayesian  simulation  could  not  be  computed  for  the 

case  when  the  shape  parameter  was  less  than  1.0  and  a 

+■  v> 
failure  existed  for  the  i   component. 

All  calculations  were  conducted  on  the  IBM  360  computer 
with  the  exception  of  the  Classic  method  with  zero  failures, 
whose  calculations  were  computed  on  a  desk  calculator. 

For  each  three-line  block  of  numbers  the  reader  may 
compare  the  three  methods'  upper  confidence  levels  com- 
puted from  the  same  arguments.   To  see  how  a  change  in  the 
scale  parameter,  the  number  of  times  a  component  fails, 
or  the  number  of  components  in  series  is  reflected  in  the 
upper  confidence  level,  the  reader  merely  moves  to  the  right 
or  down  in  the  table.   If  the  reader  wants  to  see  how  a 
change  in  the  shape  parameter  or  a  change  in  the  amount  of 
time  a  component  is  on  test  is  reflected,  he  must  go  to 
another  table. 
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V.   CONCLUSIONS 

A  "crossover"  point  in  this  discussion  will  be  defined 
as  the  point  or  general  area  at  which  an  upper  confidence 
level  that  was  previously  lower  than  one  with  which  it  was 
being  compared  becomes  higher. 

The  Bayesian  method  will  be  compared  with  the  Classic 
and  the  Semi-Classic  will  be  compared  with  Classic.   The 
Bayesian  and  the  Semi-Classic  appear  to  behave  similarly 
and  remarks  concerning  this  similarity  will  be  mentioned. 

A.   UNMODIFIED  COMPARISONS,  ALFA  =1.0 

The  first  three  tables  deal  with  a  shape  parameter  of 
1.0  and  a  system  of  components  that  either  all  of  the 
components  experience  a  failure  or  none  of  the  components 
experience  a  failure. 

In  the  cases  where  none  of  the  components  fail,  a 
crossover  point  is  exhibited  in  every  case  when  five  or 
more  components  are  in  series,  regardless  of  how  optimistic 
the  scale  parameter  becomes  or  how  long  the  test  time  is 
extended.   When  every  component  fails,  the  Bayesian  and  Semi 
Classic  systems  do  not  crossover.   In  both  cases  the  values 
of  the  upper  confidence  level  for  the  Bayesian  and  the 
Semi-Classic  methods  become  closer  together  as  the  number 
of  components  increases. 
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B.  UNMODIFIED  COMPUTATIONS,  ALFA  =0.5 

The  upper  confidence  levels  for  the  second  three  tables 
(less  the  Classic)  are  computed  with  a  more  optimistic 
prior  shape  parameter  (alfa  =  0.5)  so  naturally  the  UCL's 
for  these  two  methods  are  lower  than  those  shown  in  the 
first  three  tables. 

Crossover  points  are  much  the  same  as  in  the  case  when 
alfa  =  1.0.   The  Bayesian  and  Semi-Classic  methods  cross- 
over when  none  of  the  components  fail;  only  the  crossover 
point  occurs  when  about  10  components  are  placed  in  series. 
The  Semi-Classic  does  not  crossover  when  all  components 
fail  (the  Bayesian  is  not  examined  due  to  the  non-existance 
of  a  Gamma  random  variate  generator  with  a  non-integer 
shape  parameter) . 

C.  MODIFIED  COMPUTATIONS,  ALFA  =1.0 

The  third  three  tables  provide  a  much  more 
interesting  comparison.   Here  the  situation  where  only 
one  (or  two)  component(s)  in  the  system  fail.   The  Semi- 
Classic  and  the  Bayesian  techniques  will  crossover  the 
Classic  system  in  every  case  regardless  of  how  optimistic 
the  "priors"  are  or  how  long  the  components  are  tested. 
When  the  priors  are  more  optimistic  and/or  the  test  time  is 
larger,  then  it  takes  a  larger  number  of  components  in 
series  for  the  crossover  to  occur.   This  point  may  occur 
with  as  few  as  two  components  or  as  many  as  between  10  and 
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30  components.   The  Bayesian  UCL's  are  strictly  lower  than 
the  Semi-Classic  UCL's. 

D.   MODIFIED  COMPUTATIONS,  ALFA  =0.5 

In  the  last  three  tables  only  the  Semi-Classic  and  the 
Classic  methods  can  be  compared  for  the  Bayesian  cannot  be 
simulated  for  this  case  (non-integer). 

The  Semi-Classic  technique  is  using  a  more  optimistic 
prior  so  its  UCL's  are  lower  than  before  but  again  in  every 
case  it  will  crossover  with  the  Classic  method. 

Generally  one  could  conclude  that  as  long  as  the 
number  of  components  that  are  in  series  in  a  system  are 
few,  the  Bayesian  approach  may  yield  a  lower  value  of  an 
upper  confidence  level  for  the  system  failure  rate,  pro- 
vided the  priors  can  be  chosen  optimistically.   If  the 
system  is  complex  enough  that  more  than  "a  few"  components 
must  be  placed  in  series  then  the  Classic  system  appears  to 
yield  the  lowest  values  of  upper  confidence  levels.   The 
Semi-Classic  technique  would  only  be  useful  if  the  priors 
were  optimistic,  the  number  of  components  in  series  were 
few,  and  the  optimistic  priors  were  non-integers. 
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COMPUTER  PROGRAM  FOR  THE  CLASSIC  METHOD 


C  THE    VARIABLE    DEFINITTCNS    FOR    THIS    PROGRAM    FOLLOW: 

C  TIME  THE    NUMBER    OF    MISSION    UNITS    EACH    COMPON- 

C  ENT  IS  UNDER  TEST 

C         FAIL       THE  NUMBER  OF  TIMES  THAT  EACH  COMPONENT 

C  WILL  FAIL  DURING  THE  TOTAL  TEST  TIME 

C         XK         THE  NUMBER  OF  COMPONENTS  THAT  ARE  IN  THE 

C  SERIES  SYSTEM 

C         XLAMI      INDIVIDUAL  COMPONENT  FAILURE  RATE 

C         XLAMS      SYSTEM  FAILURE  RATE 

C         XKPRI      A  VALUE  TAKEN  FROM  TABLE  I  (BETA  VALUES) 

C 

C 

XK=1.0 

XKPRI=1. 469496 
7      TIME=30.0 
10     C=0.0 

FAIL=1.0 

XLAMUP=0.0 

XLAMI=FAIL/TIMF 

XLAMS=XK*XLAMI 

C=iXK*( XLAMI /TIME)} /XLAMS 

XLAMUP=(2.0  -XL  AMS+(XKPRI**-2)*C  +  S0RT(  (  (  4.  0*  XLAMS  )  *  (  XKPR 
1  1**2  )*C  )  +  (XKPRI**4)-*C~*2)  )/2.0 

WRITE(6,101)FAIL,XK,TIME,XLAMUP 

101  FORMAT(«  » ,5X,  • FA IL= • , F7  .3, 5X,  »XK=» ,F7.3,5X, «TIME  =  '  »F7 
2.3.5X, • XLAMUP  =  ».F10.6,/J 

TIME=TIME+20.0 

102  IF(TIMC.LE.51.0)  GO  TO  10 
TIME=TIMF+30.0 

103  IF(TIME.LE.101  .0)  GO  TO  10 
XK=XK+1.0 

XKPRI=1. 44015 

104  IF(XK.LE.2.3)  GO  TO  7 
XK=XK+2.0 

XKPRI=1. 403713 

105  IF(XK.LT.5.5)  GO  TO  7 
XK=XK+2.0 

XKPRI=1 .37725 

106  IFCXK..LT.10.5)  GO  TC  7 
XK=XK+15.0 

XKPRI=1 .34092 

107  IF( XK.LT.3D.5)  GO  TC  7 
STOP 

END 
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COMPUTER  PROGRAM  FOR  THE  SEMI-CLASSIC  METHOD 


C  THIS    SYSTEM    IS    LABELED    "SEMI-CLASSIC"    BECAUSE    IT 

C  UTILIZES    THE    "PRIORS"    THAT    WERE     INPUTS    TO    THE    BAYES- 

C  IAN    TECHNIQUE    BUT    THE    COMPUTATIONS    ARE    PERECRMED    THE 

C  THE    SAME    WAY    AS    IN    THE     "CLASSIC"    TECHNIQUE. 

C 

C      THE  VARIABLE  DEFINITIONS  EOR  THIS  PROGRAM  FOLLOW: 

C         XK         NUMBER  CF  COMPONENTS  IN  SERIES 

C         FAIL       NUMBER  CF  FAILURES  PER  COMPONENT 

C         TIME       NUMBER  OF  MISSION  UNITS  EACH  COMPONENT 

C  IS  UNCER  TEST 

C         BETA       OPTIMISTIC  SCALE  PARAMETER 

C  ALFA       OPTIMISTIC  SHAPE  PARAMETER 

C  XLAMS      SYSTEM  FAILURE  RATE 

C         XLAMI      ESTIMATE  OF  THE  COMPONENT  FAILURE  RATE 

C         XLAMUP     UPPER  CONFIDENCE  LEVEL  ON  SYSTEM 

C  FAILURE  RATE 

C 

ALFA=0.5 

6  XK=1.0 
XKPRI=1. 469496 

7  FAIL=0.0 

8  BETA=50.0 

9  TIME=30.0 

10  C=0.0 

XL AM  UP  =  0.0 

XLAMI=(ALFA+FAIL)/ (TIME+BETA) 
XLAMS=XK*XLAMI 

C=<XK*( XLAMI /{ TIME+BETA) )) /XLAMS 

XLAMUP =< (2.0^XLAMS)+( ( XKPR 1**2 )*C  )  + SORT ( ( ( 4.0*XLAMS) 
1*(  XKPRI**2KC)  +  (XKPRI»'*4)*C**2)  )/2.0 
WRITE  (  6,  10  DAL  FA,  FAIL,  XK,  TIME,  BETA,  XLAMUP 

101  FORMAT (•  • ,2Xi • ALFA=« , F4.1 ,2X, 'FAIL=' ,F4.1,2X, »XK=» 
6F5.1.2X, 'TIME='  ,F6.1,2X,'BETA=«  ,F6.1.2X,« XLAMUP=' , 
7F10.6,// ) 

TIME=TIME+20.0 

102  IF(TIME.LE.51.0)  GO  TO  10 
TIME=TIME+30.0 

103  IFCTIME.LE. 101.0)  GO  TO  10 
BETA=3ETA+50.0 

104  IF(BETA.LE. 101.0)  GO  TO  9 
FAIL=FAIL+1.0 
IF(FAIL.LE.1.1 )  GO  TO  8 
XK=XK+1.0 

XKPRI=1. 44015 

IFUK.LE.2.1  )  GO  TO  7 

XK=XK+2.0 

XKPRI=1. 403713 

IFCXK.LE.5.1 )  GO  TO  7 

XK=XK+2.0 

XKPRI=1. 37725 

IF(XK.LE.lO.l)  GO  TO  7 

XK=XK+15.0 

XKPRI=1. 34092 

IF(XK. LE.30.5)  GO  TO  7 

ALFA=ALFA+0.5 

IF1ALFA.LE.1.1 )  GO  TO  6 

STOP 

END 
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COMPUTER  PROGRAM  FOR  THE  MODIFIED  SEMI-CLASSIC  METHOD 


C  THIS    PROGRAM    IS    A    "MODIFIED    SEMI-CLASSIC"    IN    THAT 

C  THE    FAILURES    ARE    MODIFIED    SO    THAT    ONLY    ONE    COMPONENT 

C  OR    ONLY    TWC    COMPONENTS    FAIL    DURING    THE    COMPLETE    SYS- 

C  TEM    TEST.       ALL     OTHER    VARIABLES    ARE    THE    SAME    AS    FOR 

C  THE    "SEMI-CLASSIC"     PROGRAM. 

C 

C 

C  ONLY  ONE  SET  OF  VARIABLES  ARE  TESTED  WITH  THIS 

C  PROGRAM.   THE  VARIABLES  TIME , BET  A, ALFA,  AND  FAIL 

C  WOULD  HAVE  TO  BE  VARIED  TO  GENERATE  A  FULL  TABLE  OF 

C  VALUES. 

C 

c 

XKPRI=1. 44015 

XK=2.0 

TIME=100.0 

XLAMS=2. 0/100.0 

BETA=0.0 
10     ALFA=0.0 

C=1.0/TIME 

XLAMUP=O.0 

XLAMUP=(  2.0^XLAMS+(  XKPR  1**2  )-<C  + SORT (  ( ( 4. 0* XLAMS ) * ( XKPR 
1 1**2 )*C )+(XKPRI^*4)*C**2))/2.0 

WRITE( 6,101)XK,BETA ,TI ME,XLAMUP 
101    FORMATC  »,2X,«XK  =  '  ,  F7 .  3  ,  2X  ,  •  BETA  =  ',F7.3,2X, 
2«TIME  =  • ,F7.3,2X, 'XLAMUP  =  «,F10.6,/) 

XK=XK+3.0 

XLAMS=2. 0/100.0 

IFIXK.LEi5.5J  GO  TO  10 

XK=XK+2.0 

XLAMS=2. 0/100.0 

IFIXK.LE.13.5)  GC  TO  ID 

XK=XK+15.0 

XLAMS=2. 0/100.0 

IFCXK.LE.31.0)  GO  TO  10 

STOP 

END 
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COMPUTER  PROGRAM  FOR  THE  BAYESIAN  SIMULATION 
ALFA  =  1.0 


C  THE    DEFINITIONS    CF    THE    VARIABLES    USED    IN    TjHIS 

C  PROGRAM     FOLLOW: 

C 

C         TIME       THF  NUMBER  OF  MISSION  UNITS  THAT  EACH 

C  COMPONENT  IS  UNDER  TEST 

C         BETA       THE  OPTIMISTIC  SCALE  PARAMETER 

C         ALFA       THE  OPTIMISTIC  SHAPE  PARAMETER 

C         FAIL       THF  NUMBER  OF  TIMES  THAT  EACH  COMPONENT 

C  WILL  FAIL  DURING  THE  TEST  TIME 

C         K  THE  NUMBER  OF  COMPONENTS  THAT  ARE  IN  THE 

C  SERIES  SYSTEM 

C         XLAM       THE  FAILURE  RATE  FOR  AN  INDIVIDUAL 

C  COMPONENT 

C         XLAMS      THE  SYSTEM  FAILURE  RATE 

C 

C 

C      THIS  SYSTEM  IS  CALLED 

DIMENSION  XLAMS (100  0) , KEY (1000) 

TIME=30.0 

4  BETA=50.0 

5  FAIL=0.0 
7  IND  =  1 

600  CO    610     10    =     It  1000 
XLAMS(IO)     =    0.0 

610    CONTINUE 
6  ALFA=1.0 

B=TIME+BETA 

KA=ALFA+FAIL 

GO    TO    C 601, 602. 603, 60 4. 606, 607). IND 

601  K    =    1 
GO    TO    608 

602  K    =    2 
GO    TO    608 

603  K    =     5 
GO    TO    608 

6  04    K    =    10 

GO    TO    60  8 
606    K    =    30 

608  IX=999 

9  DO    50    J=l,1000 
DO    609     ILBDS    =     1,K 
TR=1.0 

10  DO    20    1=1, KA 
CALL    RANDU( IX, IY,YFL) 
IX  =  IY 

20  TR=TR*YFL 

30  XLAM=-ALGG(TR)/B 

XLAMS(J)     =    XLAMS(J)    +    XLAM 

609  CONTINUE 
KEY( J)=J 

50  CONTINUE 

CALL     SHSORT(XL AMS, KEY, 1000) 
WR ITE ( 6, 60 )K, TIME, BETA, FAIL, XLAMS (900) 
60  FORMAT!*     ',3X,'C0MP.    =     '  , I  4, 3X , • T I  ME    =    ',F7.2,3X, 

l'BETA    =     »,F7.2, 3X, • FAIL    =     ■ « F6.2. 3X. • LAMS    =    ', 
2F10.6. //) 
IND    =     IND    +    1 
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GO  TO  600 
6C7  CCNTINUE 

FAIL=FAIL+1.0 

IF(FAIL  .LE.1.5  )  GG  TO  7 

BETA=BETA+5C.O 

IF(BETA.LE. 100.1)    GO    TO    5 

TIME=TIME+2C.O 

IFCTIME.LE.50. 1)     GO    TO    4 

TIME=TTiME+30.0 

IF(TIME.LE. 100.1 )     GO   TO    4 

STOP 

END 
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COMPUTER  PROGRAM  FOR  THE 

MODIFIED  BAYESIAN  SIMULATION 

ALFA  =  1.0 


C  THIS  BAYESIAN  TECHNIQUE  IS  "MODIFIED"  IN  MUCH  THE 

C  SAME  WAY  THAT  THE  SEMI-CLASSIC  TECHNIQUE  WAS  IN  THAT 

C  THE  NUMBER  CE  FAILURES  ARE  "MODIFIED"  SO  THAT  ONLY 

C  ONE  AND  ONLY  TWO  COMPONENTS  WILL  FAIL  DURING  THE  SYS- 

C  TEM  TEST  TIME 

C 

C 

C  THE  DEFINITIONS  OF  THE  VARIABLES  FOR  THIS  PROGRAM 

C  FOLLOW: 

C         TIME       THE  NUMBER  OF  MISSION  UNITS  THAT  EACH 

C  COMPONENT  IS  UNDER  TEST 

C         BETA       THE  OPTIMISTIC  SCALE  PARAMETER 

C  K  THE  NUMBER  OF  COMPONENTS  THAT  ARE  IN 

C  THE  SERIES  SYSTEM 

C         JJ         THE  NUMBER  OF  COMPONENTS  IN  THE  SYSTEM 

C  THAT  FAIL 

C         XLAMS      THE  SYSTEM  FAILURE  RATE 

C 

c 

DIMENSION    XLAMS(IDOO) ,KEY(1000) 
1    DO    700    JJ=1,2 
TIME=30.0 

4  BFTA=50.0 

5  IND=1 

C      KA  IS  EQUAL  TO  ALFA  PLUS  FAIL,  WHICH  IS  2 

C      THE  INDEX  "J  J"  IS  EQUAL  TO  THE  NUMBER  OF  COMPONENTS 

C      THAT  FAIL  IN  THE  SERIES  SYSTEM 

600  DO  610  10=1.1000 
XLAMSI I0)=0.0 

610  CONTINUE 

B=BFTA+TIME 

GO  TO  (601, 602, 6C3, 604,606,607) ,IND 

601  K  =  l 

GO  TO  60  8 

602  K=2 

GO  TO  608 

603  K  =  5 

GO  TO  608 

604  K  =  10 

-  GO  TO  608 
606  K=30 
608  IX=999 

DO  80  J=l,1000 

DO  24  JKJ=1,JJ 

TR=1.0 
10  DO  20  I =1,2 

CALL  RANDU( I  X, I Y,YFL) 

IX=IY 
20  TR=TR*YFL 

XLAMX=-ALOG(TR)/B 

XLAMS (J )=XLAMS( J  J+XLAMX 
24  CONTINUE 

IF(K-JJ.LE.O)  GO  TO  79 

L=K-JJ 

DO  30  ILBDS=1,L 

TR=1.0 
C      KA  IS  NOW  EQUAL  TO  ONE  FOR  THE  REST  OF  THE  SYSTEM 

DO  25  I  I =1,1 

CALL  RANDU( IX, IY,YFL) 
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25 


30 
79 
80 


IX  =  IY 

TR=TR*YFL 

XLAM=-AL0G(TR)/B 

XLAMS(J)=XLAMS(J  )+XLAM 

CONTINUE 

KEYC J)=J 

CCNTINUE 

CALL  SHSORK  XL AMS , KEY , 1000 ) 

WRITE(6.90)K.TIME,BETA,XLAMS(9Q0) 


90  FORMAT(  «  '  ,3X, 
*«BETA  =  «,F7.2 
IND=INQ+1 
GC  TO  600 

607  CONTINUE 

BETA=BETA+50.0 
IF(3ETA.LE.100 
TIME=TIME+20.0 
IFCTIME.LE.50. 
TIME=TIME+30.0 
IFCTIME.LE.100 

7C0    CCNTINUE 
STOP 
END 


CCMP.    =     ' , 14, 3X, 'TIME    = 
3X,  'XLAMS    =     '  ,F10.6,//) 


1)    GO    TO    5 


•  ,F7.2,3X, 


1  )    GO    TO    4 
.1  )    GO    TO    4 


38 


COMPUTER  PROGRAM  FOR  THE  BAYESIAN  SIMULATION 
ALFA  =0.5 


c 

THIS  PROGRAM 

c 

[DISTRIBUTION 

c 

THAN  1.0 

c 

c 

c 

THE  VARIABLE 

c 

c 

K 

c 

XBETA 

c 

XTIME 

c 

c 

XFAIL 

c 

c 

NCMP 

c 

c 

XLAMS 

c 

XLAMUP 

c 

c 

c 

WILL  GENERATE  VARIATES  FROM  THE  GAMMA 
WHEN  THE  SHAPE  PARAMETER  (ALFA)  IS  LESS 


DEFINITIONS  FOR  THE  MAIN  PROGRAM  FOLLOW: 

THE  OPTIMISTIC  SHAPE  PARAMETER 

THE  OPTIMISTIC  SCALE  PARAMETER 

THE  NUMBER  OF  MISSION  UNITS  THAT  EACH 

COMPONENT  IS  UNDER  TEST 

THE  NUMBER  OF  TIMES  THAT  EACH  COMPONENT 

FAILS  CURING  TH?  TEST  TIME 

THE  NUMBER  OF  COMPONENTS  IN  THE  SERIES 

SYSTEM 

THE  SYSTEM  FAILURE  RATE 

THE  UPDER  CONFIDENCE  LEVEL  ON  THE 

FAILURE  RATE 

DIMENSION  XLAMS( 10  0  0) ,KEY(1000) 

REAL*4  K 

K  =  3.  5 

XTIME=30.0 
30  XBETA=50.0 
40  IND=1 

600  CO  610  10=1*1000 
XLAMS( I0)=0.0 

610    CONTINUE 

BETA=1 .O/CXBETA+XTIME) 

GO  T0( 6  01,602 ,6 03, 604,605* 606 ) , IND 

601  NCMP=1 

GO  TO  607 

602  NCMP=2 

GO  TO  607 

603  NCMP=5 

GO  TO  60  7  . 

604  NCMP=10 
GO  TO  607 

635  NCMP=30 
607  IX=999 

CALL  GMINIT(K,SETA) 
10  DO  50  J=l, 1000 

CO  609  ILBCS=1,NCMP 

CALL  GAMA( IX, Z) 

Z  WILL  =  AN  UNORDERED  VALUE  OF  LAMBDAI 

XLAMSU)=XLAMS  (  J  )+Z 
609    CONTINUE 

KEY( J)=J 
50    CONTINUE 

NOW  ORDER  THE  SYSTEM  FAILURE  RATES 

CALL  SHSORT( XLAMS, KEY, 1000) 

WRIT? (6, 60) NCMP, XBETA, XT  I  ME, XLAMS (900) 
60  FORMATl'  «,2X,'NCMP  = '  ,  14 , 2X , • BET A  =',F6.2,2X, 
"-•TIME  =•  ,F6.2,  2X, 'XLAMS  =',F10.6,/) 

IN0=IND+1 

GO    TO    600 
606    CCNTINUE 

XBETA=XBETA+50.0 
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IF(XRETA.LE.100. 1)     GO    TO    40 
XTIME=XTIME+20  .0 
IFCXTIME.LE.50.1)    GO    TO    30 
XTIME=XTIMc+30.0 
IFCXTIME.LE. 100.1  )     GO    TO    30 
END 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


SUBROUTINE  GMI  NI T  (  K  ,  B  ETA  ) 

ADDITIONAL  ENTRY  POINTS: 

RSULT 
GAMA(IX,Z) 


PURPOSE: 

GENERATION  OF  GAMMA  RANDOM  DEVIATES  WITH  SHAPE 
PARAMETER  LESS  THAN  ONE. 

METHOD: 

A  MODIFICATION  OF  MARSAGLIA'S  BCX-WEDGE-T A  IL  METHOD 
FOR  NORMAL  DEVIATES  IS  USED.  THE  PDF  IS  DECOMPOSED 
INTO  A  HEAD  REGION,  A  NUMBER  (DEPENDENT  ON  K)  0^ 
RECTANGLES  AND  WEDGES  AND  A  TAIL  REGION.  THE  GMINIT 
SECTION  OF  THE  SUBROUTINE  ALSO  SETS  UP  A  BINARY 
SEARCH  TREE  TO  BE  USED  FOR  EFFICIENT  SELECTION  OF  THE 
PROPER  REGION  DURING  THE  ACTUAL  GENERATING  PROCESS, 
WHICH  IS  HANDLED  BY  THE  GAMA  SECTION 

DESCRIPTION  OF  PARAMETERS 

K       GAMMA  DISTRIBUTION  SHAPE  PARAMETER 

(MUST  BE  .GE.  0.05  AND  .LE.  1.0) 
BFTA    GAMMA  DISTRIBUTION  SCALE  PARAMETER 
IX      SEED  FOR  RANDCM  NUMBER  GENERATOR 
Z       RETURNED  GAMMA  DEVIATE 

THE  PDF  OF  THE  GAMMA  FUNCTION  IS  GIVEN  BY 

F(X)  =  (1/BETA)**K  *  X**(K-1)  *  EXP (- X/BETA ) /GAMMA ( K ) 

THE  FOLLOWING  SUBROUTINES  ARE  USED: 


RANDCM( IX) 

INVGAM(K,X) 

IGAM(K.X) 


RETURNS  A  UNIFORM  (0,1)  DEVIATE 
COMPUTES  THE  INVERSE  GAMMA  CDF 
COMPUTES  THE  INCOMPLETE  GAMMA 
FUNCTICNCGAMMA  CDF) 


NOTE: 

UNDERFLOW  IS  POSSIBLE  WHEN  K  IS  LESS  THAN. 18  AND 

BECOMES  MORE  LIKELY  AS  K  DECREASES.   WHEN  K  IS  0.5 
THE  PROBABILITY  0^  UNDERFLOW  IS  ABOUT. 0C0129  FOR 
ANY  GENERATED  DEVIATE. 


C 
C 
C 
C 
C 
C 
C 
C 
C 

c 

C 
C 


SUBROUTINE    GMI  N  IT  (  K  ,  B  ET  A  ) 

REAL*4    K,INVGAM,IGAM 

INTEGER**  FIRST, TABLE, BOTTOM, END 

LOGICAL^l  USED 

DIMENSION  P( 10 0),X( 101), H(100), 0(100), R(100),  9(100) 

CI  MENS  ION  TABLE(202),PR0B(202),NEXT(202) , LAST (202) 

DIMENSION  TEST (202) .LIST (202 ), USE D( 20 2) 

DIMENSION  RAND(2) 

EQUIVALENCE  (U,RAND(1)),(V,RAND(2)) 

THIS  FIRST  SECTION  INITIALIZES  CONSTANTS  AND  TABLES  TO 

BE  USED  BY  GAMA  WHEN  IT  IS  CALLED.   THE  FOLLOWING 
ARE  USED  BY  GAMA: 

PO  PROBABILITY  FOR  "HEAD"  REGION 

PN  PROBABILITY  FOR  "TAIL"  REGION 

P(I)  PROBABILITY  FCR  I-TH  RECTANGLE 

X(I)  LEFT-HAND  BOUNDARY  OF  I-TH  RECTANGLE 

H(I)  WIDTH  OF  I-^H  RECTANGLE 

0(1)  PROBABILITY  CF  I-TH  WEDGE 
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c 

R<  I) 

c 

B(  I) 

c 

ALPHA 

t 

PROS 

c 

TABLE 

c 

c 

FIRST 

c 

NEXT, 

c 

LAST 

c 

Jl 

c 

HI  TO 

c 

H4 

c 

c 

c 

CHECK 

REJECTION  TEST  RATIO  FOR  I-TH  WEDGE 
Y  INTERCEPT  FOR  I-TH  WEDGE 
SHAPE  PARAMETER  -  1. 
ORDERED  VECTOR  CF  PROBABILITIES 

VECTOR  OF  WEDGE/RECTANGLE  NUMBERS  CORRESPOND- 
ING TO  PRCBAEILITIESS  IN  PROB 
STARTING  POINT  FOR  BINARY  SEARCH 
LINKS  FOR  BINARY  SEARCH 

POSITION  IN  PROB  OF  PO 
CONSTANTS  FOR  APPROXIMATION  TO  INVERSE  GAMMA 
CDF  FOR  SMALL  VALUES  OF  Z 

CHECK  FOR  K  IN  RANGE 

IF( (K.GE.0.05)  .AND.  { K .LE. 1 .0 ) ) GO  TO  5 

WRITE{6,4)K 

4  FCRMAT( //' OGMI NIT  CALLED  WITH  K=',1PE16.6, 
*•   CUT  OF  RANGE'/) 

RETURN 
C 
C      GET  UPPER  BOUND  ON  NUMBER  OF  RECTANGLES 

5  N=20.+6.6/K 
IF(N.GT.100)N=100 
N=2*N+2 

NM=M-1 
ALPHA=K-1. 
GK=GAMMA(KJ 
P3=5.E-5/(K*K) 
HFAC=2. 
C 

C      SET  UP  RECTANGLE  BOUNDS 
C 

X(1)=1NVGAM(K, PO) 
P3=IGAM(K,X(1) ) 
H(1)=.25*X(1 ) 
CO  10  I=2«N 

X( I )=X( I-l)+H( 1-1) 
H(  I  )  =  H( T-l )*HFAC 
P(I )=0. 
0( I )=0. 
10  CONTINUE 

X(N+1) =X(N)+H(N) 
C 

C      ZERO  PROBABILITY  VECTORS  AND  LINKS 
DO  15  I =  1,M 

NEXT( I )=0 
LAST( I )=0 
PROB (I ) =0. 
LISTd  )=0 
USED( I )=. FALSE. 
15  CONTINUE 
C 

C      FIND  COEFFICIENTS  FOR  NEWTON-R APHSON  APPROXIMATION  TO 
C      SLOPE  OF  PDF 
C 

B1=-2.*ALPHA 
B2=ALPHA*(ALPHA-1.  ) 
A1=B1+1. 

A2=ALPHA*( ALPHA-2.  ) 
C=l. -ALPHA 
C 

C      FIND  RECTANGLE  PROBABILITIES  AND  WEDGE  VALUES 
C 

PL=PO 

FL=EXP(ALPHA*AL0G(X(1) )-X( 1) )/GK 

DO  40  I =  1.N 

FU=EXP(ALPH ANALOG (X( I+1))-X(I+1))/GK 
PU=IGAM(K,X( 1+1) ) 
P( I )=H( I)*EU 
0(  I  )  =  PU-PL-P(I  ) 
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C  NEWTON-RAPHSON  ITERATION  TO  FIND  POINT  WHERE 

C  SLOPE  OF  PDF  IS  PARALLEL  TO  CHORD 

C 

W  =  X(I ) 

S=(FU-FL)  /H(I  ) 

SC=S*GK 

DO  20  J  =  l,8 

Y=W*( ( W+A1)*W+A2+SC*EXP(C*AL0G(W)+W) )/( (W+Bl 
IFIASS(Y-W) .LT.1.E-4*Y)G0  TO  30 
W  =  Y 
20       CONTINUE 
C 

C  FIND  VALUES  FCR  REJECTION  METHOD 

C 

30       B( I  )=EXP( ALPHA*ALOG(Y)-Y)/GK+S*(X(  I)-Y) 
R(I )=(B(I )-FU)/(FL-FU) 
C 
C  TEST  TO  SEE  IF  ENOUGH  RECTANGLES  HAVE  BEEN  TAKEN 

IF(PU.GT. 0.999)G0  TO  45 
C 

C  RESET    PROBABILITY    AND    FUNCTION    VALUES    FOR    NEXT 

C  RECTANGLE 

C 

PL  =  PU 
FL=FU 
40    CONTINUE 
C 
C  FIND    LOWER    END    OF    NCN-ZERO    PROBABILITY    VECTOR 

45    L0W=2*( N-I )+l 
C 

C  FIND    TAIL    PROBABILITY 

SUM1=P0 
DO    50    J=1,N 

SUM1=SUM1+P( J)+0( J) 
50    CONTINUE 
PN=1.-SUM1 
C 

C  GFT    VALUES    TO    COMPUTE     LOW    END    APPROXIMATION    TO    INVERSE 

C  GAMMA    CDF 

C 

H1=K*GK*P0 
H2=l./K 

H3=(K+1. )*5.0E-7 
H4=4./(K+1. ) 
C 

C      GENERATE  PROBABILITY  VECTOR 
C 

DO  80  1=1, N 

PROB(  I  )=P ( I  ) 
TABLEU  )  =  I 
PROB( I  +  N) =0(1  ) 
TABLE( I+N )=-I 
80    CONTINUE 

PROB(M-l)=P0 
TABLE(M-1)=0 
PROB(M)=PN 
TABLE(M) =0 
C 

C  SORT    PROBAEILITY    VECTOR 

C 

DO    110    1=1, MM 
ICH=0 
L  =  M-I 
DO    100    J=1,L 

IF(PROB(J)-PROB( J+l) ) 100 » 100,90 
90  TEMP=PROB(J) 

PROBt J)=P3C8 (J+l) 
PR0B(J  +  1  )=TFMP 
ITEMP=TAELE( J) 
TABLE( J) =TABLE( J+l) 
TABLE( J+1)=ITEMP 
ICH=1 
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c 
c 
c 
c 


c 
c 
c 
c 


100 
110 


120 


130 


140 


150 


1  51 


152 
153 


154 
155 


156 
157 


1571 

158 
1525 


159 


1591 
160 

161 


162 
163 


CONTINUE 

IF(  ICH) 120, 120,110 
CCNTINUE 
PROB (M)  =  l. 

CCNVFRT  PRCB  TO  CUMULATIVE  PROBABILITIES  AND 
FIND  FIRST  AND  Jl 

J1  =  0 

FIRST=0 
L=L0W+1 
DO  130  I=L.M 

IF( (TABLE( I ) .FO.O)  .AND.  ( PROB ( I ) . E 0. PO ) ) Jl  =  1 

PR0BCI)=PROB(  D+PR0BCI-1J 

IF( (PRCB( I) .GE.0.5)     .AND.     { F I RST . EQ.O  )  ) F  I  RST= I 
CONTINUE 
IF(FIRST.EO.M)FIRST=MM 

IF(J1)14C,140, 150 
J1=L0W 


NOW  DETERMINE 
SEARCH 


THE  VECTORS  NEXT  AND  LAST  FOR  BINARY 


BOTTOM= 
END  =  1 
PR=.25 
LIST(l) 
TFSTd  ) 
USED(FI 
DO  159 
LI 
IF 
PR 
DO 

CO 
IF 
J  = 

IF 
NF 
LI 
IF 
PR 
TF 
DO 

CO 
IF 
J  = 

IF 
LA 

GO 
J  = 

EN 
LI 
TE 
CCNTINU 
BOTTCM= 
PR  =  PR*0 
1  =  1 

IF{LIST 
BOTTCM= 
IF(BOTT 
DC  162 
LI 
TF 
CONTINU 
GO  TO  1 
USEDtLI 
1=1  +  1 
IFU.LE 


=  FI 
=  .5 
RST 
1=1 
=  LI 
(US 
N=T 
15 

NTI 
(  .N 
J-l 
(LI 
XT( 
ST( 
(  (L 
L=T 
ST( 
15 

NTI 
(  .N 
J  +  l 
(LI 
ST( 
TO 
LI 
D  =  E 
ST( 
ST( 

END 
.  5 


RST 

)=.TRUE. 

.BOTTOM 

ST(I  ) 

EC(L 1  +  1) )G0  TO  155 

EST(  I  )+PR 

2  J = LOU, MM 

IF(PROB( J) .GT.PRNJGO  TO  153 

NUC 

OT.  USED( J) )GC  TO  154 

-J)153,155,155 

LI  )  =  J 

I  )=NIEXT(LI) 

I.EO.LOW)  .OR.  USED(LI-l)  )  GO  TO  159 

EST(I)-PR 

I  )=PRN 

6  J=LOW,MN 

IF(PROB(  J)  .GT.PRDGO  TO  157 

NUE 

OT.  USED( J) )GC  TO  1571 

-JJ158, 158,157 

LI  )  =  J 
1585 

ND+1 
END)=J 
END) =PRL 


( I) )160,160. 163 

BOTTOM- 1 

OM) 165,165, 161 

J  =  1  ,  BOTTOM 

ST( J) =LIST( J+l) 

ST( J)=TEST( J+l) 

c 

591 

ST(  I  )  )=.TRUE. 

.BOTTOMJGO  TO  1591 
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END=BOTTCM 

GO  TO  151 
C 

C      SETUP  ^IRST  CALL  TO  RANDOM 

165  CALL  OVFLOW 
C 

C      THROUGH  WITH  SETUP  PHASE.  QUIT. 

RETURN 
C 

c 

C      THIS  SECTIGN  PRINTS  OUT  THE  VALUES  GENERATED  BY  GMINIT 
C 

ENTRY  RSULT 

WRITF(6,166)K,  SETA 

166  FORMATt/'IGENERATED  VALUES  FOR  K=',1PE14.6, 
*•    BETA  =  «  ,E14.6) 

WRITE(6, 170) P0,PN 
170  FORMAT! «  OHE AD  PROB ABI L I TY= • , 1PE 1 5. 6 , 
*■    TAIL  PR03A3ILITY=' , E15.6) 

WRITE (6, 18  0) 
180  FORMAT!  ' OR ECTANGL E /WEDGE  VALUE S ' //2X , '  I  •  , 9X , ' X ( I ) «  , 12X 
*  , «  P(I )«  ,  12X,'Q< I  )'  ,12X, «R( I) ',12X,  'B( I) '/) 

SUM1=0. 

SUM2=0. 

CO  192  1=1, N 

IF(P( I  )  )  193,193,185 
185  SUM1=SUM1+P( I) 

SUM2=SUM2+0(I) 

WRITP(6,190)I,X(I),P{I),0(I),R(I),B(I) 
190  FORMAT (1X13, 1P5E16. 6) 

192  CONTINUE 

193  WRITE(6,194) SUM1.SUM2 

194  FORMAT* 'OSUMS' ,15X,1D2E16.6) 
WRITE(6,195)J1,H1,H2,H3,H4 

195  FORMAT  (/ 'OVAL'JES  FOR  HEAD/TAIL  APPROXIMATION:'/'  Jl=«, 
*I6/'  H1  =  ',E16.6,'  F2  =  ',E16.6,«  H3=',E16.6,'  H4«, 
1E16.6) 

WRITE(6,196)FIRST 

196  FORMAT(/'OSTART!NG    POINT    FOR    BINARY    SEARCH', 14) 
WRITE (6,  197) (I  , PROB(I)  , TA&L E ( I  )  ,NE XT( I )  ,  L A  ST ( I )  ,  I =1 , M) 

197  FORMAT (/'0PR03ABILITY  VECTOR  AND  L I NKS ' //  14X,  » PRGB » , 9X 
*,  'TABLE'  ,4X,'NEXT',5X,'LAST'//(1XI5,1PE16.6,0P3I9)} 

RETURN 
C 

C      THIS  IS  THE  SECTION  THAT  ACTUALLY  COMPUTES  THE  RANDOM 
C      VARIATES 
C 

ENTRY  GAMAHX,  Z) 
C 
C      GET  TWO  UNIFORM  DEVIATES 

CALL  RANDOM <IX,U,2) 
C 
C      CONDUCT  BINARY  SEARCH  OF  PROBABILITY  VECTOR 


C 


c 


J  =  F 


200  IF(U-PROB( J) )210,25C,230 
210  IF(L4ST( J) )250,250,220 
220  J=LAST(J) 
GO  TO  200 
230  IF(NEXT{ J) ) 245 ,245,240 
2  40  J=NEXT{ J ) 
GO  TO  200 
245  J=J+1 
C 

C      LOCATED  PROBABILITY  DIVISION.  GET  TABLE  VALUE. 
250  N=TABLE(J) 

IF(N)260,290,320 
C 
C      THIS  SECTION  IS  FOR  THE  WEDGES 


260  N=-N 

CALL  RANDOM! IX,U,1) 
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c 
c 
c 


270  IF(U.LE.V)GC  TO  280 
TFMP=U 

u=v 

V=TEMP 
280  Z=X(N)+H(N)*U 

IF( V.LE.R(N) )G0  TO  330 

THIS  STEP  IS  PERFORMED  VERY  RARELY 
W=U+EXP(£LPHA*ALOG(Z)-Z)/B(N) 
IFCV.LE.WJGO  TO  33  0 
CALL  RAMD0VUX.U.2) 
GO  TO  270 


290 


C 
C 

c 

c 
c 
c 


THIS  SECTION  IS  FOR  HEAD/TAIL  PROBABILITIES 


IFU.EQ.J1  JGO  TO  300 
Z=INVGAM(K,1.-PN*V) 
GO  TO  330 
300  Z=(H1*V)**H2 

IF(Z.LT.H3)G0  TO  330 

Z  =  2.*Z/(  l.  +  SORT(  l.-H4*Z) ) 

GO  TO  330 

THIS  SECTION  IS  FOR  THE  RECTANGLES 

320  Z=X{N)+H(N)*V 

SCALE  GAMMA  VARIATE  AND  EXIT 

330  Z=Z*BETA 
RETURN 
END 
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INVGAM    SUBROUTINE 


FUNCTION     INVGAM(K,Z) 

RFAL*4    INVGAM,IGAM,K,EPS/l.E-08/ 

IF(Z.GT.0.0)GO    TO    10 

INVGAM=0.0 

RFTURN 

10   x=z 

t=K-l. 

G=GAMMA(K) 

DO    40    1=1,30 

Y  =  G*( IGAM(K,X)-Z)/ (  EXP  ( T* ALOG(  X  )-X)  ) 
20    IF(Y.LT.X)GC    TO    30 

Y=Y*.5 

GO    TO    20 
30    X=X-Y 

IF(ABS(Y)-EPS*X)50,5C,40 
40    CONTINUE 
50    INVGAM=X 

RFTURN 

THIS  IS  THE  END  OF  THE  INVGAM  SUBROUTINE 

END 
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IGAM  SUBROUTINE 


10 

20 

40 


FUNCTION  IGAM(K.X) 

IMPLICIT  REAL*8  (D) 

REAL*4  ICAM.K 

REAL*8  EPS/1. 0-13/ 

IFCX.GT.O. )GQ  TO  5 

IGAM=0. 

RETURN 

IFfX.LE. 12.0)G0  TO  8 

IGAM=1.0 

RETURN 

DX=D3LF( X) 

CK=DBLE(K) 

DTERM=DX**DK/DGAMMA (CK) 

DSUM=DTERM/DK 

CO  20  1=1.30 

IF (DABS ( DTERM) -EPS*DSUM) 40,4 0,10 

DTERM=DTERM*{-OX )/DcLOAT( I ) 

CSUM=DSUM+DTERM/ (DK+CFLOAT ( I ) ) 

CONTINUE 

IGA,M  =  SNGL(  DSUM  ) 

R-TURN 

THIS  IS  THE  END  01=  THE  IGAM  SUBROUTINE 

END 
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